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ABSTRACT 

We study planetesimal evolution in circumbinary disks, focusing on the three systems Kepler 16, 34 
and 35 where planets have been discovered recently. We show that for circumbinary planetesimals, 
in addition to secular forcing, eccentricities evolve on a dynamical timescale, which leads to orbital 
crossings even in the presence of gas drag. This makes the current locations of the circumbinary 
Kepler planets hostile to planetesimal accretion. We then present results from simulations including 
planetesimal formation and dust accretion, and show that even in the most favourable case of 100% 
efficient dust accretion, in situ growth starting from planetesimals smaller than ~ 10 km is difficult 
for Kepler 16b, Kepler 34b and Kepler 35b. These planets were likely assembled further out in the 
disk, and migrated inward to their current location. 
Subject headings: planets and satellites: formation 



1. INTRODUCTION 

In the last decade planets have been found in very per- 
turbed systems such as close binary star systems. The 
first of these p lanets to be discovered were orbiting the 



primary star (Queloz et al. 2000 Hatzes et al.| 2003 



Zucker et al.|20Tl4[ ), but the latest additions to the fa mily, 



after promi sing results using stellar eclipse timings (Lee 

KI 7 " 



et al. 



2009) , involve plane ts in circumbinary orbits: 



pleOFtDoyle et al.||2011[ ) and Kepler 34 and 35 QWelsh 



et al. 2012p . The parameters of these new planets are 



summarised in Table [T] 

The existence of planets in these systems baffles planet 
formation theory. A crucial step in the process of build- 
ing a planet, namely growing gravitationally bound pro- 
toplanets from km-sized planetesimals, can be hindered 
or stopped in these perturbed en vironments for planetes- 
imals on circumprimary orbits ([Marzari fc Scholl||2000[ 
Thebault et al.|2006[ [Paardekooper et al.|2008||Thebault| 



201ip . The coupling between gravitational perturba- 
tions of the companion star and gas drag stirs up the 
eccentricities of planetesimals, which leads to high en- 
counter velocities. This makes accretion towards larger 
bodies difficult. Similar problems haunt planetesimals on 



circumbinary orbit sdMoriwaki fc Nak agawa 2004 ; |Scholl 
et al.|2007| |jviarzari et al.||2008| |Mes chiari 20T3T~^ 



Above studies focused on gravitational dynamics and 
gas drag only. In this work, we investiga te the effect of 
collisions on t he ev olution of the system. Paardekooper] 
& Leinhardt (2010) showed that in a system with high- 
speed collisions, it is necessary to keep track of collision 
outcomes. Notably, if collisions are mostly destructive, 
any surviving planetesimals are embedded in a sea of 
small debris. If they can pick up some of this debris, 
planetesimals can grow despite th e hostile environ ment 
dPaardekooper & Leinhardt||2010l |Xie et al.||2010h. In 



1 DAMTP, Wilberforce Road, Cambridge CB3 OWA, United 
Kingdom; S.Paardekooper@damtp.cam.ac.uk 

2 School of Physics, University of Bristol, H.H. Wills Physics 
Laboratory, Tyndall Avenue, Bristol BS8 1TL, U.K. 

3 Observatoire de Paris, F-92195 Meudon Principal Cedex, 
France 



TABLE 1 

Binary and planet parameters 

Kepler 16 a Kepler 34 6 Kepler 35 b 



M A /M G 


0.69 


1.0 


0.89 


M B /M @ 


0.20 


1.0 


0.81 


a b /AU 


0.22 


0.23 


0.18 


eb 


0.16 


0.52 


0.14 


Mp/Mj 
a p /o b 


0.33 


0.22 


0.13 


3.2 


4.7 


3.4 


e p 


0.0069 


0.18 


0.042 



a [Doyle et al.| pOll]) 
b [Welsh et al.| { |2012) 

this Letter, we aim to explore this possibility in the newly 
found planet-harbouring systems of Kepler 16, 34 and 35. 

We begin in Section [2] by reviewing the eccentricity 
evolution of planetesimals in circumbinary orbits. We 
discuss the model in Section [3| present the results in 
Section [4] and conclude in Section [5] 

2. SECULAR AND NON-SECULAR ECCENTRICITY 
EVOLUTION 

Consider the dynamics of massless particles at semi- 
major axis a and period 2n/fl around a binary with 
masses Ma and Mb, total mass M* = Ma + Mb, or- 
bital period 27r/f2h, eccentricity eb and semi-major axis 
Ob- Secular perturbation theory gives an evolution equa- 
tion for the com plex eccentricity E of planete simals in 
the gas-free case (Moriwaki & Nakagawa 1 |2004| ): 

1 d 2 E 

where r is the secular timescale and Ef — ef exp(izuh) is 
the complex forcing, with z&b the longitude of periastron 
of the binary orbit and 



e f 
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3el/4 



Seg/2 



(2) 



the forced eccentricity ( |Moriwaki fc Na kagawa 2004h . 
The secular timescale is given by (Moriwaki & Nak agawa] 



2 
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Equation ([I]) describes an oscillation around the forced 
eccentricity with an amplitude given by \E(t = 0) — J5f|. 
The period of oscillation is the secular timescale, which is 
longer than the dynamical timescale for a 3> &b- While 
initially, planetesimal orbits are phased, the spatial fre- 
quency of the oscillations increases with time, so that 
eventually orbital crossings occur lead ing to high en- 
counter velocities ( Thebault et al.|2006 ). In the presence 
of a gas disk and associated drag forces, these oscilla- 
tions are damped, and size-dep endent equilibrium orbits 
exist ( |Paardekooper et ah 2008 ) . Even if orbital crossings 
can be prevented, the size-dependence of equilibrium or- 
bits leads to high encounter velocities between bodies of 
different size. This is called differential orbital phasing 
(Thebault et al. 2006). For equal-mass binaries, such as 
Kepler 34 and 35, ef ~ 0. While this is favourable for 
accretion, it is not the whole story. 

Secular perturbation theory is valid on timescales 
longer than both the binary period and the local orbital 
timescale. In addition, planetesimal eccentricities evolve 
on a shorter timescale, e ven in the case eb = 0, wher e 
secular effects are absent (Moriwaki & Nakagawa 2004). 
This short timescale evolution can be obtained by aver- 
aging the disturbing function over the mean longitude of 
the binary orbit only, and expanding terms up to second 
order in both eb and planetesimal eccentricity e. The 
eccentricity can then be shown to oscillate, on a local 
orbital timescale, around an eccentricity 



e ff = 7 



3 M A M B (a^z 
4 M? 



U 



1 



34 „ 



(4) 



where the subscript ff indicates a fast timescale. Note 
that, unlike secular oscillations, the period is indepen- 
dent of eb, and that eg falls off faster with distance than 
ef. In the case of almost equal mass binaries (like Ke- 
pler 34 and 35), ef will be very small, while eff can be 
significant. Since e now evolves on a local dynamical 
timescale, orbital crossings will occur after only a few 
local orbits. Gas drag is unable to damp these fast oscil- 
lations, because it acts on longer timescales than a dy- 
namical timescale for km-sized objects and realistic gas 
densities. Orbit crossings cannot be prevented. Typical 
encounter velocities will be ~ egaQ, (250 m/s in Kepler 
16, and 500 m/s in Kepler 34 and 35 at the locations 
of the planets) . Despite the lack of secular forcing, this 
makes the systems Kepler 34 and 35 hostile to accre- 
tion at the current planet positions. For Kepler 16, en- 
counter velocities around 250 m/s can lead to accretion 
only w hen planetesimals hav e reached sizes of 50 km or 
larger (Thebault et al.|2006). At twice the current semi- 
T6E7 



major axis of Kepler 16b, encounter velocities go down a 
factor of 4, making accretion possible at that location if 
pl anetesimals can somehow grow to > 10 km (also seen 
by |Meschiari]|2012[ ). 



3. MODEL INGREDIENTS 

In our simulations we consider a system of two stars 
with a coplanar circumbinary disk. The gas component 



of the disk is assumed to be circular and orbiting the 
binary centre of mass. The solid component of the disk 
consists of planetesimals > 1 km, which we model as 
particles, and small dust, on the same orbits as the gas. 
Planetesimals can form from small dust, accrete small 
dust on their surface, and be returned to dust in catas- 
trophic collisions. Below, we explain the planetesimal 
evolution model in more detail. 

3.1. Two-dimensional approximation 

As in Paardekooper & Leinhardt (2010), we restrict 
ourselves to planetesimal orbits lying in the orbital plane 
of the binary. Preliminary calculations allowing for incli- 
nation of planetesimal orbits show no qualitative differ- 
ences. In addition, allowing for three-dimensional (3D) 
motions requires an impractical increase in CPU time. 
Moreover, we expect collision velocities induced by fast 
or secular forced eccentricities to be much higher than 
the escape velocity from the surface of the planetesi- 
mal (catastrophic collisions), which, as we explain be- 
low, makes the motions essentially two-dimensional (2D). 
However, confining all orbits to a single plane will under- 



Leinhardt 


2010) 


system, w 


lich i 



evolution of the 



namical and collisional effects. In a strongly perturbed 
system, where almost every collision is destructive, it is 
difficult to estimate the 3D collision timescale since it is 
unclear what the inclination distribution of the fragments 
will be. F or catastrophic collisions, Leinhardt fc Stew-] 



art] (2012) showed that the maximum velocities obtained 
by the largest fragments are comparable to the escape 
speed of the combined projectile-target mass. Therefore, 
in the regime of high-speed collisions expected in a close 
binary system the planetesimal disk is expected to stay 
approximately as thin as in the unperturbed case (as- 
suming the largest collision remnants trace the majority 
of the mass). Thus, the thickness of the disk is set by 
the escape velocity of the planetesimals. In 2D, we can 
vary the collision timescale by changing the total mass in 
solids (results do not depend sensitively on this mass). 
Note that the thin disk approximation does not neces- 
sarily apply to the small debris. If the dust disk thickens 
because of collisions both dust accretion and planetesi- 
mal formation will be effected. 

3.2. Planetesimal formation and orbital phasing 

Planetesimal dynamics in binary systems is strongly 
affected by the timescale on which planetesimals are 
formed. If planetesimals form fast compared to the 
timescale of eccentricity forcing, and in a single burst, the 
interaction with the binary stars will keep their orbits in 
phase and collision velocities low ( Heppenheimer||l978 ). 
However, size-dependent gas drag introduces differential 
orbital phasing, which leads to high collision velocities 
between bodies of different size, maki ng planetesimal ac- 
cretion difficult ( Thebault et al.|2"0 06). Eccentricity forc- 
ing towards eff oc^urs~61i^^y^ianiicaI timescale. On such 
short timescales orbital crossings can not be prevented by 
gas drag. In addition, if planetesimals form continuously, 
their orbits will be out of phase from the start. The local 
collision velocity is estimated as, v co i ~ eaQ, where e is 
a typical eccentricity. 



3 



A reasonable timescale for the formation o f an individ- 
ual planetesimal is 10 4 years 
a lower limit on the timesca 



Lissauer 1993 1 . This places 
e tor dust to be converted 



into planetesimals. Chambers (2010) found that the lat- 
ter timescale could vary by several orders of magnitude 
depending on the local conditions in the disk. In the 
simulations presented here, as in |Paardekooper fc Lein- 



hardt (2010), planetesimals form continuously with half 
of the total (local) dust mass converted into planetesi- 
mals in 10 5 local orbits. Planetesimals that wander off 
the computational domain are added to the innermost 
or outermost dust bin, thereby allowing them to be re- 
cycled. 

3.3. Collisions 
In previous, purely dynamical codes, planetesimals 



drag (e.g. 


Marzari & 


Scholl et a 


. 2007}-^ 



■fe Zhou||2009| |Meschiari||2012p In the work presented" 
here this is not necessarily the case, because collisions 
can affect the size of planetesimals on their way to their 
steady-state orbits. 

If we assume that the thickness of the planetesimal 
disk is not strongly affected by high-speed collisions (Sec. 
3.1) and a typical collision speed v co \ = eafi, we can 
estimate the collision timescale in a 3D disk. Consider 
a population of planetesimals of size R, mass M, and 
number density n. The collision timescale is: 

n 



0.5 / a \ : 
T VAU/ 



R V 17 g/cm 2 
km J £ 



Pp 



3 g/cm 3 



irR?nv coi 

'Me 



(5) 



where we have used n = £/(2AzM), with S the surface 
density of planetesimals and a thickness Az set by the 
escape velocity, Az = v csc /(2a£l). A similar expression 
can be derived in the 2D approximation, yielding a ratio: 



Z^E =2 073 (-^) 



3/2 



(6) 



which depends only on a. To get realistic collision 
timescales in our 2D approximation, we can tune the T2D 
value to T3D by artificially reducing the disk mass by a 
factor 2000. This setting is, to first order, an acceptable 
one because we have found that the qualitative outcome 
of the simulations does not depend sensitively on the disk 
mass. 

Having chosen a disk mass, we want to represent this 
mass as accurately as possible. We cannot track all indi- 
vidual bodies of km-size since their number can exceed 
10 1 1 . Instead, we inflated the radius of each particle (see 
e.g. Thebault fc Brahi"c||1999 ). We found that taking an 
inflated radius for a 1 km planetesimal, R lri i/a < 10~ 4 , 
gave satisfying collision statistics without introducing a 
bias in encounter velocities. Collision outcomes are de- 
termined based on t he velocity-dependent catastr ophic 
disruption criteria of |Stewart fc Leinhardt ( |2009[ ) (see 
Paardekooper & Leinhardt 20 1T^ 
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3.4. Gas and dust disk 
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Fig. 1. — Planetesimal eccentricity and semi-major axis for bi- 
nary parameters of Kepler 16AB. Top panel: no dust accretion, 
t = 250,000 binary orbits. Bottom panel: full dust accretion, 
t = 80, 000 binary orbits. Colour indicates the size of the planetes- 
imal in km. All darker-coloured particles have accreted mass since 
the start of the simulation. The vertical dashed line indicates the 
position of Kepler 16b. The vertical dotted l ine the inner boundary 
of the accretion-friendly zone identified by |Marzari et al.| l |2012[ ). 
The solid lines indicate 2eg (steep) and 2ef. 

The gas disk is assumed to be static with no pres- 
sure gradient so that there is no radial drift for solids 
that are coupled to the gas through aerodynamic drag. 
While radial drift could be easily incorporated for the 
planetesimals it will be more important for the small 
dust component, especially for mass that is in m-sized 
bodies. Since we have no information about the size dis- 
tribution of bodies < 1 km, we choose to neglect radial 
drift. For simplicity, the dust surface density is oc r , 
so that the mass inside a ring of radial thickness Ar 
is constant throughout the disk. Experiments with dif- 
ferent density profiles showed no qualitative differences 
in the outcome of the simulations. Radial migration of 
planetesimals therefore does not play a major role in our 
simulations. This may not be true for the small dust 
component, whose radial drift will be strongly affected 
by pressure structure in the disk. 

3.5. Dust accretion 

If most collisions between planetesimals are destruc- 
tive, a large fraction of the solid component of the disk 
will be in small debris, which can be picked up by remain- 
ing planet esimals. This process involves an e fficiency fac- 
tor ed (see Paardekooper & Leinhardt 2010). In this let- 
ter, we consider the two extreme cases ed = (no dust 
accretion ) and en = 1 (full dust accret ion). In view of the 
results of Housen & Holsapple ( 2011[ ), we switch off dust 
accretion (even when = 1) whenever the relative ve- 
locity between planetesimal and dust exceeds 100 times 
the escape velocity from the surface of the planetesimal. 

4. RESULTS 

We focus first on the Kepler 16 system. We take a com- 
putational domain for the planetesimal disk 2.5 < r/a D < 
50, with a total solid mass of 1.5 • 10 8 Mi, where Mi is 
the mass of a 1 km planetesimal, and a surface density 
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Fig. 2. — Histogram of collision velocities for Kepler 16 scenario 
without dust accretion, t = 250, 000 binary orbits. Colours indicate 
different radii, with r in units of a^. 




Fig. 3. — Mass distribution in units of 1 km planetesimal for 
the simulations depicted in Fig. [T] Dotted lines: without dust 
accretion. Solid lines: dust accretion. Green indicates planetesimal 
mass, red dust mass, and black total mass. Mass is binned radially 
into 32 bins of equal size, so that initially, M(r) is constant because 
£ <x r _1 . 
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Initially, all the mass is in small dust. Planetesi 



mals form at a r ate given by e p = 10 ( Paardekooper & 



Leinhardt 2010[ ) and have an initial size ol 1 km. The gas 
disk is assumed to be circular, with a (constant) density 
of 1.4 • 10~ 9 g/cm 3 . 

4.1. No dust accretion 

First we consider a case without dust accretion. 
Growth can then only occur through collisions. In the 
top panel of Fig. [l]we show the eccentricities of planetes- 



imal orbits versus semi-major axis. For the binary pa- 
rameters of Kepler 16, ef dominates over eg throughout 
the disk. While gas drag acts to slowly force planetes- 
imals onto equilibrium orbits, the continuous introduc- 
tion of new planetesimals with e = and the collisional 
destruction of older planetesimals lead to orbits that are 
unphased. Collisions between planetesimals on unphased 
orbits would occur at v co i = efafl w 21 (at, /a) 3 / 2 km/s. 
This is in good agreement with what is measured from 
the simulations, see Fig. [2j We find that accretion is 
possible for a/a^ > 20, or a > 4.4 A U, i n good agree 



ment with results of Meschiari| (2012) and Marzari et al. 
(2012|. Even though many collisions around a/a^ = 20 
are still destructive, enough accreting collisions occur to 
sustain a population of larger planetesimals. 

4.2. Full dust accretion 

We now turn on dust accretion in full, which is the 
most favourable case possible for accretion. Dust accre- 
tion is quenched for relative velocities exceeding 100i> CS c- 
The bottom panel of Fig. [T] shows that in this case it 
is possible to grow planetesimals further in. However, 
even in this most favourable case we find it impossible 
to grow planetesimals at the location where Kepler 16b 
resides (vertical dashed line in Fig. [IJ. In regions where 
accretion is possible, a/ab > 6, all the mass ends up in 
planetesimals, and the maximum planetesimal size that 
can be reached is only limited by the amount of mass 
available. This is illustrated in Fig. [3j which shows that 
while for ed = less than 1% of the total mass ends up 
in planetesimals, for ed — 1 essentially all the mass in- 
side r/ab = 20 resides in planetesimals. After 250,000 
binary orbits, the region r/a^ > 20 has not evolved to 
the planetesimal-only state yet. Figure [3] also shows an- 
other important effect, namely mass transport towards 
accretion- friendly regions. The region inside r/a^ = 6 is 
depleted, not only in planetesimals because of destruc- 
tive collisions, but also of small dust. This can happen 
because planetesimals are not necessarily formed and de- 
stroyed at the same radius. If a planetesimal is born at 
radius r with e = and its eccentricity is excited to e = e 
its orbit lies between r(l — e) and r(l + e). The plan- 
etesimal can be destroyed at any radius between these 
extremes. At the inner edge of the accretion-friendly 
region, mass gets locked up in large, indestructible plan- 
etesimals. Outward mass transport into the accretion- 
friendly region is therefore a one-way street, which even- 
tually leads to the depletion of the inner disk. This pro- 
cess competes against inward radial dust drift, which we 
neglect here. 

4.3. Other systems 

We now briefly discuss the systems Kepler 34 and Ke- 
pler 35. These are both almost equal-mass binaries (see 
Table JT]), which makes ef small. Planetesimal orbital 
evolution is predominantly due to short timescale inter- 
actions with the binary, whose amplitude falls off rapidly 
with distance, since eg oc a~ 2 . This makes these systems 
slightly more accretion friendly than Kepler 16 at small 
radii. In the case of no dust accretion, we find that plan- 
etesimal growth is possible for a/a^ > 12 for Kepler 34 
and for a/a^ > 15 for Kepler 35. Full dust accretion 
leads to very similar results: planetesimal accretion is 
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possible slightly further in compared to Kepler 16. How- 
ever, in situ accretion starting from 1 km planetesimals 
is still not possible. 

5. DISCUSSION AND CONCLUSIONS 

We have studied planetesimal collisions in circumbi- 
nary gas disks, focusing on the planet-harbouring sys- 
tems Kepler 16, 34 and 35. We have shown that in ad- 
dition to secular forcing, planetesimals experience eccen- 
tricity forcing on a dynamical timescale, which leads to 
eccentricity oscillations and orbital crossings that can not 
be prevented by gas drag. This makes the current loca- 
tion of the planets Kepler 16b, 34b and 35b very hostile 
for planetesimal accretion. 

We then used a numerical model similar to that of 



Paardekooper & Leinhardt ( 2010[ ) including planetesimal 
formation and accretion of small dust. Even in the most 
favourable case of 100% efficient dust accretion, we have 
been unable to grow planetesimals from initially 1 km at 
the current location of the planets. Since dust accretion 
is likely to be less than 100% efficient, for example be- 
cause not all the small dust will be concentrated in the 
midplane of the disk, we conclude that in situ planetesi- 
mal accretion is difficult for the planets Kepler 16b, 34b 
and 35b. 

We have made several necessary simplifications to 
make following the collisional evolution of the planetes- 
imal population tractable. First of all, we have ig- 
nored gas dynamics throughout and worked with a static 
circular gas disk. While the gas disk is likely to be- 
co me eccentric, especially at small radii, it was shown 



Paardekooper et al. ( 2008 ) that unless the gas relaxes 



towards the forced eccentricity, including gas dynamics 
makes matters worse for planetesimal accretion. For the 
fast eccentricity oscillations to be damped by gas drag, 
the gas disk will have to oscillate in phase with the plan- 
etesimals. Full hydrodynamical simulations are neces- 
sary to determine whether this is the case. These can 
also be used to study the effect of the inner truncation 



of the gas disk, and we will consider such simulations in 
future investigations. 

We considered the planar case, but a small inclination 
of the binary plane with respect to the gas disk may pro- 
mote planetesimal accretion in the circumprimary case 
( Xie fc Zhou||2009| . However, because of the fast eccen- 
tricity oscillations and the resulting orbital crossings it is 
unclear if this effect can play a ro le in the circumbinary 
case. Moreover, it was shown in Fragner et al. (20111 
that including gas dynamics again makes matters worse, 
even in the misaligned disk case. 

A formation mechanism which can leapfrog the prob- 
lematic km-size range, such as gravitational collaps e 
aided by streaming instabilities (Johansen et al. 2007), 
may overcome the problems of planetesimal accretion. It 
remains to be seen, however, if such a mechanism can op- 
erate in close binary systems. Preliminary calculations 
show that in the current model, we would need to start 
with planetesimals of at least 10 km in order for in situ 
accretion of the Kepler circumbinary planets to become 
possible. 

The most straightforward solution is that the three 
circumbinary planets were assembled further out in an 
accretion-friendly region, and migrated in towards their 
current location at a later stage. This can be achieved at 
a relatively early stage, in the 10-100 km size range, by 
radial drift due to a pressure gradient in the gas, or at a 
later stage when the planet is more or less fully grown, 
by Type I or Type II planetary migration. Whatever the 
migration mechanism, it is likely that the inner edge of 
the truncated gas disk will cause migration to stall. We 
then expect the current location of the planets to be close 
to the truncation radius of the gas disk. 



We thank the referee, H. Perets, for an insightful re- 
port. SJP, ZML and CB acknowledge support from 
STFC. 
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